Open In Colab

Dyadic Representational Similarity Analysis (DyRSA)

The following code implements the method described in:

Rhoads, S. A. & Marsh, A. A. (2023). Modeling the minds of friends: Dorsomedial prefrontal cortex encodes accurate dyadic meta-evaluations about how others perceive us. PsyArXiv. https://doi.org/10.31234/osf.io/g9e3z
In [1]:
%matplotlib inline
%reload_ext autoreload
%autoreload 2
In [ ]:
import numpy as np, pickle
import matplotlib.pyplot as plt, seaborn as sns
from scipy.stats import spearmanr
from urllib.request import urlopen

import nibabel as nib
from nilearn.image import new_img_like
from nilearn.plotting import view_img_on_surf
from nltools.stats import fdr

from joblib import Parallel, delayed
from tqdm import tqdm
from copy import deepcopy

import dyrsa
In [6]:
# Load participant information
partner_dict = {'sub-1001':'sub-2001','sub-2001':'sub-1001',
                'sub-1002':'sub-2002','sub-2002':'sub-1002',
                'sub-1003':'sub-2003','sub-2003':'sub-1003',
                'sub-1004':'sub-2004','sub-2004':'sub-1004',
                'sub-1005':'sub-2005','sub-2005':'sub-1005',
                'sub-1006':'sub-2006','sub-2006':'sub-1006',
                'sub-1007':'sub-2007','sub-2007':'sub-1007',
                'sub-1008':'sub-2008','sub-2008':'sub-1008',
                'sub-1009':'sub-2009','sub-2009':'sub-1009',
                'sub-1010':'sub-2010','sub-2010':'sub-1010',
                'sub-1011':'sub-2011','sub-2011':'sub-1011',
                'sub-1012':'sub-2012','sub-2012':'sub-1012',
                'sub-1013':'sub-2013','sub-2013':'sub-1013',
                'sub-1015':'sub-2015','sub-2015':'sub-1015'}

dyad_ids = [1, 1, 2, 2, 3, 3, 4, 4, 
            5, 5, 6, 6, 7, 7, 8, 8, 
            9, 9, 10, 10, 11, 11, 
            12, 12, 13, 13, 14, 14]

subject_ids = ['sub-1001','sub-2001',
               'sub-1002','sub-2002',
               'sub-1003','sub-2003',
               'sub-1004','sub-2004',
               'sub-1005','sub-2005', 
               'sub-1006','sub-2006', 
               'sub-1007','sub-2007', 
               'sub-1008','sub-2008', 
               'sub-1009','sub-2009', 
               'sub-1010','sub-2010', 
               'sub-1011','sub-2011', 
               'sub-1012','sub-2012', 
               'sub-1013','sub-2013', 
               'sub-1015','sub-2015']

nsubj = len(subject_ids)
In [7]:
# create friend matrix and nonfriend_matrix
friend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
friend_matrix[:] = np.nan
nonfriend_matrix = np.zeros((len(subject_ids), len(subject_ids)))
nonfriend_matrix[:] = np.nan

for i, sub_i in enumerate(subject_ids):
    for j, sub_j in enumerate(subject_ids):
        if (sub_i == partner_dict[sub_j]) or (partner_dict[sub_i] == sub_j):
            friend_matrix[i, j] = 1
            friend_matrix[j, i] = 1
        else:
            nonfriend_matrix[i, j] = 1
            nonfriend_matrix[j, i] = 1
        
        # also exclude self
        if i == j:
            nonfriend_matrix[i, j] = np.nan
            friend_matrix[i, j] = np.nan

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 6))

sns.heatmap(friend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[0], cbar=False, cmap='RdBu')

axs[0].set_title('Friend Pairs')

sns.heatmap(nonfriend_matrix, xticklabels=subject_ids, yticklabels=subject_ids, ax=axs[1], cbar=False, cmap='RdBu')
axs[1].set_title('Non-friend Pairs')

plt.show()
In [8]:
# Load dictionary of dictionary of searchlight objects contains RDMs based on the cross-validated Euclidean distance between the neural patterns across trials and runs
# `dyad_rdms` is formatted like: dyad_rdms[sub_id][cond_id]
# e.g., {'sub-1001' : {'Meta': searchlight_object, 'Other': searchlight_object}, 
#        'sub-1002' : {'Meta': searchlight_object, 'Other': searchlight_object}, 
#        ...}
# pickled file is available at https://osf.io/t5cmd/
try:
    with open('../data/dyad_rdms.pickle', 'rb') as handle:
        dyad_rdms = pickle.load(handle)
except:
    # if not downloaded, load from OSF
    dyad_rdms = pickle.load(urlopen("https://osf.io/t5cmd/data/dyad_rdms.pickle", 'rb'))
In [9]:
# Load mask
mask_fn = 'MNI152NLin2009c-Schaefer2018-DMN21-k100_3mm.nii.gz'
mask_file = nib.load(mask_fn)
mask_img = mask_file.get_fdata()
mask = mask_img==1
x, y, z = mask_img.shape
In [ ]:
# Run DyRSA for every possible pair, then can filter friend vs non-friend pairs later
# first grab nvoxels from first subject
voxel_indices = list(dyad_rdms['sub-1001']['Meta'].rdm_descriptors['voxel_index'])
nvoxels = len(voxel_indices)

# initialize result matrices
pairwise_zmat_brain = np.empty((nvoxels, nsubj, nsubj)) 
pairwise_zmat_brain[:] = np.nan

# loop through all pairs and compute DyRSA
for subject_index, subject_id in tqdm(enumerate(subject_ids)):        
    for subject_index2, subject_id2 in enumerate(subject_ids):
        if subject_index == subject_index2:
            pairwise_zmat_brain[:, subject_index, subject_index2] = np.nan
        else:
            # Neural Meta(A(B(A))) ~ Neural Other(B(A))
            subj_sl_RDM = dyad_rdms[subject_id]['Meta']
            partner_sl_RDM = dyad_rdms[subject_id2]['Other'] #can also replace with other model RDMs
            slresults = Parallel(n_jobs=-1)(delayed(spearmanr)(partner_sl_RDM.dissimilarities[x_index], 
                            x.dissimilarities[0,]) for x_index, x in enumerate(subj_sl_RDM))
            spearmanrhos_brain = [float(e[0]) for e in slresults]
            pairwise_zmat_brain[:, subject_index, subject_index2] = np.array([np.arctanh(r) for r in deepcopy(spearmanrhos_brain)])
In [20]:
# remove friend values and flatten across voxels to shape nvox, sum(nonfriend_matrix)=728
all_nonfriend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0], 
                                                  int(np.sum(nonfriend_matrix==1))))

# remove non-friend values and flatten across voxels to shape nvox, nsubj
all_friend_pairs_flat = np.zeros((pairwise_zmat_brain.shape[0], nsubj))

# loop through voxels
for v in range(pairwise_zmat_brain.shape[0]):
    all_nonfriend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][nonfriend_matrix==1]
    all_friend_pairs_flat[v] = pairwise_zmat_brain[v,:,:][friend_matrix==1]
In [ ]:
# compute the mean of all_friend_pairs_flat
actual_group_mean = np.nanmean(all_friend_pairs_flat, axis=1) # mean at every voxel

# generate null distribution of bootstrap samples
nboot = 5000
zstat, perm_p = dyrsa.run_permutations(nsubj, nboot, 
                                        actual_group_mean, 
                                        all_nonfriend_pairs_flat,
                                        replacement=True,
                                        write_zstat=True)
In [ ]:
# Write unthresholded result
results_brain = np.zeros((x*y*z))
results_brain[voxel_indices] = zstat
results_brain = results_brain.reshape((x, y, z))
unthresh_img = new_img_like(mask_file, results_brain)
In [45]:
# plot
view_img_on_surf(unthresh_img,threshold=1)

Out[45]:
In [ ]:
# Apply FDR correction and write thresholded result
mask_pval = [pv for mask_ind, pv in zip(mask.ravel(), perm_p) if mask_ind]
fdr_thresh = fdr(np.array(mask_pval), q=.05)

results_brain2 = np.zeros((x*y*z))
results_brain2[voxel_indices] = np.array([m if p<fdr_thresh else np.nan for m, p in zip(zstat, perm_p)])
results_brain2 = results_brain2.reshape((x, y, z))
thresh_img = new_img_like(mask_file, results_brain2)
In [43]:
# plot
view_img_on_surf(thresh_img,threshold=.001)
Out[43]:
In [ ]: